Make the prediction grids
Max Lindmark
07 March, 2023
Intro
Make an evenly spaced UTM prediction grid with all spatially varying covariates for the diet and the biomass data
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(sp)
library(raster)
library(devtools)
library(RCurl)
library(sdmTMB)
library(terra)
library(ncdf4)
library(chron)
# Source code for map plots
source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5
# Source code for lon lat to utm
source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/lon-lat-utm.R")
theme_set(theme_plot())Read data and depth-raster
# Read data
d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") %>%
rename(X = x, Y = y)
#> Rows: 9376 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (3): country, haul_id, ices_rect
#> dbl (20): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed 2 variables (X, Y)Make the grid with depth
First make a grid for the biomass data, then subset that based on the extend of the stomach data
x <- d$X
y <- d$Y
z <- chull(x, y)
coords <- cbind(x[z], y[z])
coords <- rbind(coords, coords[1, ])
plot(coords[, 1] ~ coords[, 2]) # plot data
sp_poly <- sp::SpatialPolygons(
list(sp::Polygons(list(sp::Polygon(coords)), ID = 1))
)
sp_poly_df <- sp::SpatialPolygonsDataFrame(sp_poly,
data = data.frame(ID = 1)
)
cell_width <- 3
pred_grid <- expand.grid(
X = seq(min(d$X), max(d$X), cell_width),
Y = seq(min(d$Y), max(d$Y), cell_width),
year = unique(d$year)
)
ggplot(pred_grid %>% filter(year == 2019), aes(X, Y)) +
geom_point(size = 0.1) +
theme_void() +
coord_sf()
#> filter: removed 850,905 rows (96%), 31,515 rows remaining
sp::coordinates(pred_grid) <- c("X", "Y")
inside <- !is.na(sp::over(pred_grid, as(sp_poly_df, "SpatialPolygons")))
pred_grid <- pred_grid[inside, ]
pred_grid <- as.data.frame(pred_grid)
ggplot(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000)) +
geom_point(size = 0.001, alpha = 0.5) +
NULL
#> filter: removed 512,784 rows (96%), 18,992 rows remaining
plot_map +
geom_point(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000), size = 0.001, alpha = 0.5) +
NULL
#> filter: removed 512,784 rows (96%), 18,992 rows remaining
#> Warning: Removed 704 rows containing missing values (geom_point).
# Add lat and lon
# Need to go from UTM to lat long for this one...
# https://stackoverflow.com/questions/30018098/how-to-convert-utm-coordinates-to-lat-and-long-in-r
xy <- as.matrix(pred_grid %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000))
#> mutate: changed 531,776 values (100%) of 'X' (0 new NA)
#> changed 531,776 values (100%) of 'Y' (0 new NA)
v <- vect(xy, crs="+proj=utm +zone=33 +datum=WGS84 +units=m")
y <- project(v, "+proj=longlat +datum=WGS84")
lonlat <- geom(y)[, c("x", "y")]
pred_grid$lon <- lonlat[, 1]
pred_grid$lat <- lonlat[, 2]
ggplot(filter(pred_grid, year == 1999), aes(lon, lat)) + geom_point()
#> filter: removed 512,784 rows (96%), 18,992 rows remaining
# Add depth now to remove islands and remaining land
# https://gis.stackexchange.com/questions/411261/read-multiple-layers-raster-from-ncdf-file-using-terra-package
# https://emodnet.ec.europa.eu/geoviewer/
dep_raster <- terra::rast("data/Mean depth natural colour (with land).nc")
class(dep_raster)
#> [1] "SpatRaster"
#> attr(,"package")
#> [1] "terra"
crs(dep_raster, proj = TRUE)
#> [1] "+proj=longlat +datum=WGS84 +no_defs"
plot(dep_raster)
pred_grid$depth <- terra::extract(dep_raster, pred_grid %>% dplyr::select(lon, lat))$elevation
ggplot(pred_grid, aes(lon, lat, color = depth*-1)) +
geom_point()
pred_grid$depth <- pred_grid$depth*-1
pred_grid <- pred_grid %>% drop_na(depth)
#> drop_na: removed 74,340 rows (14%), 457,436 rows remaining
pred_grid %>%
filter(year == 1999) %>%
drop_na(depth) %>%
#mutate(water = ifelse(depth < 0.00000001, "N", "Y")) %>%
ggplot(aes(X*1000, Y*1000, fill = depth)) +
geom_raster() +
NULL
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> drop_na: no rows removed
plot_map +
geom_point(data = pred_grid, aes(X*1000, Y*1000), size = 0.001) +
geom_sf()
#> Warning: Removed 19628 rows containing missing values (geom_point).
plot_map +
geom_raster(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000, fill = depth), size = 0.001) +
geom_sf()
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> Warning: Ignoring unknown parameters: size
#> Warning: Removed 762 rows containing missing values (geom_raster).
Substrate
First add lat and lon based on X and Y (utm)
# # Read raster
# substrate <- raster("data/substrate_tif/BALANCE_SEABED_SEDIMENT.tif")
#
# # Reproject the raster to fit the UTM pred grid...
# sr <- "+proj=utm +zone=33 +datum=WGS84 +units=m "
#
# # Project Raster... This takes some time
# projected_sub_raster <- projectRaster(substrate, crs = sr)
#
# plot(projected_sub_raster)
#
# # Now extract the values from the saduria raster to pred_grid
# pred_grid$substrate <- raster::extract(projected_sub_raster, pred_grid %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000))
#
# unique(pred_grid$substrate)
#
# # Plot
# ggplot(pred_grid, aes(X, Y, color = substrate)) +
# geom_point()
#
# factor(sort(unique(round(pred_grid$substrate))))
#
# pred_grid$substrate <- round(pred_grid$substrate)
#
# pred_grid <- pred_grid %>% mutate(substrate = ifelse(substrate == 1, "bedrock", substrate),
# substrate = ifelse(substrate == 2, "hard-bottom complex", substrate),
# substrate = ifelse(substrate == 3, "sand", substrate),
# substrate = ifelse(substrate == 4, "hard clay", substrate),
# substrate = ifelse(substrate == 5, "mud", substrate))
#
# pred_grid <- pred_grid %>% drop_na(substrate)
#
# # I. Bedrock.
# # II. Hard bottom complex, includes patchy hard surfaces and coarse sand (sometimes also clay) to boulders.
# # III. Sand including fine to coarse sand (with gravel exposures).
# # IV. Hard clay sometimes/often/possibly exposed or covered with a thin layer of
# # sand/gravel.
# # V. Mud including gyttja-clay to gyttja-silt.
#
# # Plot
# ggplot(pred_grid, aes(X, Y, fill = substrate)) +
# geom_raster() +
# coord_sf()Oxygen
For these variables, we add quarter to the prediction grid
# Add quarter
pred_grid_1 <- pred_grid %>% mutate(quarter = 1)
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA
pred_grid_4 <- pred_grid %>% mutate(quarter = 4)
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA
dat <- bind_rows(pred_grid_1, pred_grid_4)
# Downloaded from here: https://resources.marine.copernicus.eu/?option=com_csw&view=details&product_id=BALTICSEA_REANALYSIS_BIO_003_012
# Extract raster points: https://gisday.wordpress.com/2014/03/24/extract-raster-values-from-points-using-r/comment-page-1/
# https://rpubs.com/boyerag/297592
# https://pjbartlein.github.io/REarthSysSci/netCDF.html#get-a-variable
# Open the netCDF file
ncin <- nc_open("data/NEMO_Nordic_SCOBI/dataset-reanalysis-scobi-monthlymeans_1664182224542.nc")
print(ncin)
#> File data/NEMO_Nordic_SCOBI/dataset-reanalysis-scobi-monthlymeans_1664182224542.nc (NC_FORMAT_CLASSIC):
#>
#> 1 variables (excluding dimension variables):
#> float o2b[longitude,latitude,time]
#> long_name: Sea_floor_Dissolved_Oxygen_Concentration
#> missing_value: NaN
#> standard_name: mole_concentration_of_dissolved_molecular_oxygen_in_sea_water
#> units: mmol m-3
#> _FillValue: NaN
#> _ChunkSizes: 1
#> _ChunkSizes: 523
#> _ChunkSizes: 383
#>
#> 3 dimensions:
#> time Size:336
#> axis: T
#> long_name: Validity time
#> standard_name: time
#> units: days since 1950-01-01 00:00:00
#> calendar: gregorian
#> _ChunkSizes: 512
#> _CoordinateAxisType: Time
#> valid_min: 15721.5
#> valid_max: 25917.5
#> latitude Size:523
#> axis: Y
#> standard_name: latitude
#> long_name: latitude
#> units: degrees_north
#> _CoordinateAxisType: Lat
#> valid_min: 48.49169921875
#> valid_max: 65.8914184570312
#> longitude Size:383
#> standard_name: longitude
#> long_name: longitude
#> units: degrees_east
#> axis: X
#> _CoordinateAxisType: Lon
#> valid_min: 9.01375484466553
#> valid_max: 30.2357654571533
#>
#> 24 global attributes:
#> references: http://www.smhi.se
#> institution: Swedish Meterological and Hydrological Institute
#> history: See source and creation_date attributees
#> Conventions: CF-1.5
#> contact: servicedesk_cmems@mercator-ocean.eu
#> comment: Provided by SMHI as a Copernicus Marine Environment Monitoring Service production unit
#> bullentin_type: reanalysis
#> cmems_product_id: BALTICSEA_REANALYSIS_BIO_003_012
#> title: CMEMS V4 Reanalysis: SCOBI model 3D fields (monthly means)
#> FROM_ORIGINAL_FILE__easternmost_longitude: 30.2357654571533
#> FROM_ORIGINAL_FILE__northernmost_latitude: 65.8914184570312
#> FROM_ORIGINAL_FILE__westernmost_longitude: 9.01375484466553
#> FROM_ORIGINAL_FILE__southernmost_latitude: 48.49169921875
#> shallowest_depth: 1.50136542320251
#> deepest_depth: 711.059204101562
#> source: SMHI reanalysis run NORDIC-NS2_1d_20201201_20201201
#> file_quality_index: 1
#> creation_date: 2021-11-09 UTC
#> bullentin_date: 20201201
#> start_date: 2020-12-01 UTC
#> stop_date: 2020-12-01 UTC
#> start_time: 00:00 UTC
#> stop_time: 00:00 UTC
#> _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
# Get longitude and latitude
lon <- ncvar_get(ncin,"longitude")
nlon <- dim(lon)
head(lon)
#> [1] 9.013755 9.069310 9.124865 9.180420 9.235975 9.291530
lat <- ncvar_get(ncin,"latitude")
nlat <- dim(lat)
head(lat)
#> [1] 48.49170 48.52503 48.55836 48.59170 48.62503 48.65836
# Get time
time <- ncvar_get(ncin,"time")
time
#> [1] 15721.5 15751.0 15780.5 15811.0 15841.5 15872.0 15902.5 15933.5 15964.0
#> [10] 15994.5 16025.0 16055.5 16086.5 16116.0 16145.5 16176.0 16206.5 16237.0
#> [19] 16267.5 16298.5 16329.0 16359.5 16390.0 16420.5 16451.5 16481.0 16510.5
#> [28] 16541.0 16571.5 16602.0 16632.5 16663.5 16694.0 16724.5 16755.0 16785.5
#> [37] 16816.5 16846.5 16876.5 16907.0 16937.5 16968.0 16998.5 17029.5 17060.0
#> [46] 17090.5 17121.0 17151.5 17182.5 17212.0 17241.5 17272.0 17302.5 17333.0
#> [55] 17363.5 17394.5 17425.0 17455.5 17486.0 17516.5 17547.5 17577.0 17606.5
#> [64] 17637.0 17667.5 17698.0 17728.5 17759.5 17790.0 17820.5 17851.0 17881.5
#> [73] 17912.5 17942.0 17971.5 18002.0 18032.5 18063.0 18093.5 18124.5 18155.0
#> [82] 18185.5 18216.0 18246.5 18277.5 18307.5 18337.5 18368.0 18398.5 18429.0
#> [91] 18459.5 18490.5 18521.0 18551.5 18582.0 18612.5 18643.5 18673.0 18702.5
#> [100] 18733.0 18763.5 18794.0 18824.5 18855.5 18886.0 18916.5 18947.0 18977.5
#> [109] 19008.5 19038.0 19067.5 19098.0 19128.5 19159.0 19189.5 19220.5 19251.0
#> [118] 19281.5 19312.0 19342.5 19373.5 19403.0 19432.5 19463.0 19493.5 19524.0
#> [127] 19554.5 19585.5 19616.0 19646.5 19677.0 19707.5 19738.5 19768.5 19798.5
#> [136] 19829.0 19859.5 19890.0 19920.5 19951.5 19982.0 20012.5 20043.0 20073.5
#> [145] 20104.5 20134.0 20163.5 20194.0 20224.5 20255.0 20285.5 20316.5 20347.0
#> [154] 20377.5 20408.0 20438.5 20469.5 20499.0 20528.5 20559.0 20589.5 20620.0
#> [163] 20650.5 20681.5 20712.0 20742.5 20773.0 20803.5 20834.5 20864.0 20893.5
#> [172] 20924.0 20954.5 20985.0 21015.5 21046.5 21077.0 21107.5 21138.0 21168.5
#> [181] 21199.5 21229.5 21259.5 21290.0 21320.5 21351.0 21381.5 21412.5 21443.0
#> [190] 21473.5 21504.0 21534.5 21565.5 21595.0 21624.5 21655.0 21685.5 21716.0
#> [199] 21746.5 21777.5 21808.0 21838.5 21869.0 21899.5 21930.5 21960.0 21989.5
#> [208] 22020.0 22050.5 22081.0 22111.5 22142.5 22173.0 22203.5 22234.0 22264.5
#> [217] 22295.5 22325.0 22354.5 22385.0 22415.5 22446.0 22476.5 22507.5 22538.0
#> [226] 22568.5 22599.0 22629.5 22660.5 22690.5 22720.5 22751.0 22781.5 22812.0
#> [235] 22842.5 22873.5 22904.0 22934.5 22965.0 22995.5 23026.5 23056.0 23085.5
#> [244] 23116.0 23146.5 23177.0 23207.5 23238.5 23269.0 23299.5 23330.0 23360.5
#> [253] 23391.5 23421.0 23450.5 23481.0 23511.5 23542.0 23572.5 23603.5 23634.0
#> [262] 23664.5 23695.0 23725.5 23756.5 23786.0 23815.5 23846.0 23876.5 23907.0
#> [271] 23937.5 23968.5 23999.0 24029.5 24060.0 24090.5 24121.5 24151.5 24181.5
#> [280] 24212.0 24242.5 24273.0 24303.5 24334.5 24365.0 24395.5 24426.0 24456.5
#> [289] 24487.5 24517.0 24546.5 24577.0 24607.5 24638.0 24668.5 24699.5 24730.0
#> [298] 24760.5 24791.0 24821.5 24852.5 24882.0 24911.5 24942.0 24972.5 25003.0
#> [307] 25033.5 25064.5 25095.0 25125.5 25156.0 25186.5 25217.5 25247.0 25276.5
#> [316] 25307.0 25337.5 25368.0 25398.5 25429.5 25460.0 25490.5 25521.0 25551.5
#> [325] 25582.5 25612.5 25642.5 25673.0 25703.5 25734.0 25764.5 25795.5 25826.0
#> [334] 25856.5 25887.0 25917.5
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
#> [1] 336
tunits
#> $hasatt
#> [1] TRUE
#>
#> $value
#> [1] "days since 1950-01-01 00:00:00"
# Get oxygen
dname <- "o2b"
oxy_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(oxy_array)
#> [1] 383 523 336
# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")
# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])
# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))
# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)
# Replace netCDF fill values with NA's
oxy_array[oxy_array == fillvalue$value] <- NA
# Next, we need to work with the months that correspond to the quarters that we use.
# loop through each time step, and if it is a good month save it as a raster.
# First get the index of months that correspond to Q4
months
#> [1] 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1
#> [26] 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2
#> [51] 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3
#> [76] 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4
#> [101] 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5
#> [126] 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6
#> [151] 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7
#> [176] 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8
#> [201] 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9
#> [226] 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10
#> [251] 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11
#> [276] 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12
#> [301] 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1
#> [326] 2 3 4 5 6 7 8 9 10 11 12
index_keep_q1 <- which(months < 4)
index_keep_q4 <- which(months > 9)
oxy_q1 <- oxy_array[, , index_keep_q1]
oxy_q4 <- oxy_array[, , index_keep_q4]
months_keep_q1 <- months[index_keep_q1]
months_keep_q4 <- months[index_keep_q4]
years_keep_q1 <- years[index_keep_q1]
years_keep_q4 <- years[index_keep_q4]
# Now we have an array with data for that quarter
# We need to now calculate the average within a year.
# Get a sequence that takes every third value between 1: number of months (length)
loop_seq_q1 <- seq(1, dim(oxy_q1)[3], by = 3)
loop_seq_q4 <- seq(1, dim(oxy_q4)[3], by = 3)
# Create objects that will hold data
dlist_q1 <- list()
dlist_q4 <- list()
oxy_1 <- c()
oxy_2 <- c()
oxy_3 <- c()
oxy_ave_q1 <- c()
oxy_10 <- c()
oxy_11 <- c()
oxy_12 <- c()
oxy_ave_q4 <- c()
# Now average by quarter. The vector loop_seq_q1 is 1, 4, 7 etc. So first i is 1, 2, 3,
# which is the index we want.
for(i in loop_seq_q1) { # We can use q1 as looping index, doesn't matter!
oxy_1 <- oxy_q1[, , (i)]
oxy_2 <- oxy_q1[, , (i + 1)]
oxy_3 <- oxy_q1[, , (i + 2)]
oxy_10 <- oxy_q4[, , (i)]
oxy_11 <- oxy_q4[, , (i + 1)]
oxy_12 <- oxy_q4[, , (i + 2)]
oxy_ave_q1 <- (oxy_1 + oxy_2 + oxy_3) / 3
oxy_ave_q4 <- (oxy_10 + oxy_11 + oxy_12) / 3
list_pos_q1 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)
list_pos_q4 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)
dlist_q1[[list_pos_q1]] <- oxy_ave_q1
dlist_q4[[list_pos_q4]] <- oxy_ave_q4
}
# Now name the lists with the year:
names(dlist_q1) <- unique(years_keep_q1)
names(dlist_q4) <- unique(years_keep_q4)
# Now I need to make a loop where I extract the raster value for each year...
# Filter years in the pred_grid based on quarter
d_sub_oxy_q1 <- dat %>% filter(quarter == 1)
#> filter: removed 457,436 rows (50%), 457,436 rows remaining
d_sub_oxy_q4 <- dat %>% filter(quarter == 4)
#> filter: removed 457,436 rows (50%), 457,436 rows remaining
# Create data holding object
oxy_data_list_q1 <- list()
oxy_data_list_q4 <- list()
# Create factor year for indexing the list in the loop
d_sub_oxy_q1$year_f <- as.factor(d_sub_oxy_q1$year)
d_sub_oxy_q4$year_f <- as.factor(d_sub_oxy_q4$year)
# Loop through each year and extract raster values for the data points
for(i in sort(unique(d_sub_oxy_q1$year_f))) { # We can use q1 as looping index, doesn't matter!
# Set plot limits
ymin = 54; ymax = 59; xmin = 12; xmax = 22
# Subset a year
oxy_slice_q1 <- dlist_q1[[i]]
oxy_slice_q4 <- dlist_q4[[i]]
# Create raster for that year (i)
r_q1 <- raster(t(oxy_slice_q1), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
r_q4 <- raster(t(oxy_slice_q4), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
# Flip...
r_q1 <- flip(r_q1, direction = 'y')
r_q4 <- flip(r_q4, direction = 'y')
plot(r_q1, main = paste(i, "Q1"))
plot(r_q4, main = paste(i, "Q4"))
# Change projection to UTM (same as pred grid)
sr <- "+proj=utm +zone=33 +datum=WGS84 +units=m "
proj_raster_q1 <- projectRaster(r_q1, crs = sr)
proj_raster_q4 <- projectRaster(r_q4, crs = sr)
# Filter the same years and select only coordinates
d_slice_q1 <- d_sub_oxy_q1 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)
d_slice_q4 <- d_sub_oxy_q4 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)
# Make into a SpatialPoints object
data_sp_q1 <- SpatialPoints(d_slice_q1, proj4string = CRS(sr))
data_sp_q4 <- SpatialPoints(d_slice_q4, proj4string = CRS(sr))
# Extract raster value (oxygen)
rasValue_q1 <- raster::extract(proj_raster_q1, data_sp_q1, method = "bilinear")
rasValue_q4 <- raster::extract(proj_raster_q4, data_sp_q4, method = "bilinear")
# Add in the raster value in the df holding the coordinates for the data
d_slice_q1$oxy <- rasValue_q1
d_slice_q4$oxy <- rasValue_q4
# Add in which year
d_slice_q1$year <- i
d_slice_q4$year <- i
# Now the unit of oxygen is mmol/m3. I want it to be ml/L. The original model is in unit ml/L
# and it's been converted by the data host. Since it was converted without accounting for
# pressure or temperature, I can simply use the following conversion factor:
# 1 ml/l = 103/22.391 = 44.661 μmol/l -> 1 ml/l = 0.044661 mmol/l = 44.661 mmol/m^3 -> 0.0223909 ml/l = 1mmol/m^3
# https://ocean.ices.dk/tools/unitconversion.aspx
d_slice_q1$oxy <- d_slice_q1$oxy * 0.0223909
d_slice_q4$oxy <- d_slice_q4$oxy * 0.0223909
# Create a index for the data last where we store all years (because our loop index
# i is not continuous, we can't use it directly)
index_q1 <- as.numeric(as.character(d_slice_q1$year))[1] - 1992
index_q4 <- as.numeric(as.character(d_slice_q4$year))[1] - 1992
# Add each years' data in the list
oxy_data_list_q1[[index_q1]] <- d_slice_q1
oxy_data_list_q4[[index_q4]] <- d_slice_q4
}
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
# Now create a data frame from the list of all annual values
big_dat_oxy_q1 <- dplyr::bind_rows(oxy_data_list_q1)
big_dat_oxy_q4 <- dplyr::bind_rows(oxy_data_list_q4)
big_dat_oxy <- bind_rows(mutate(big_dat_oxy_q1, quarter = 1),
mutate(big_dat_oxy_q4, quarter = 4))
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA
#> mutate: new variable 'quarter' (double) with one unique value and 0% NATemperature
# Open the netCDF file
ncin <- nc_open("data/NEMO_Nordic_SCOBI/dataset-reanalysis-nemo-monthlymeans_1664183191233.nc")
print(ncin)
#> File data/NEMO_Nordic_SCOBI/dataset-reanalysis-nemo-monthlymeans_1664183191233.nc (NC_FORMAT_CLASSIC):
#>
#> 1 variables (excluding dimension variables):
#> float bottomT[longitude,latitude,time]
#> standard_name: sea_water_potential_temperature_at_sea_floor
#> units: degrees_C
#> long_name: Sea floor potential temperature
#> missing_value: NaN
#> _FillValue: NaN
#> _ChunkSizes: 1
#> _ChunkSizes: 523
#> _ChunkSizes: 383
#>
#> 3 dimensions:
#> time Size:336
#> axis: T
#> long_name: Validity time
#> standard_name: time
#> units: days since 1950-01-01 00:00:00
#> calendar: gregorian
#> _ChunkSizes: 512
#> _CoordinateAxisType: Time
#> valid_min: 15721.5
#> valid_max: 25917.5
#> latitude Size:523
#> axis: Y
#> standard_name: latitude
#> long_name: latitude
#> units: degrees_north
#> _CoordinateAxisType: Lat
#> valid_min: 48.49169921875
#> valid_max: 65.8914184570312
#> longitude Size:383
#> standard_name: longitude
#> long_name: longitude
#> units: degrees_east
#> axis: X
#> _CoordinateAxisType: Lon
#> valid_min: 9.01375484466553
#> valid_max: 30.2357654571533
#>
#> 24 global attributes:
#> references: http://www.smhi.se
#> institution: Swedish Meterological and Hydrological Institute
#> history: See source and creation_date attributees
#> Conventions: CF-1.5
#> contact: servicedesk_cmems@mercator-ocean.eu
#> comment: Provided by SMHI as a Copernicus Marine Environment Monitoring Service production unit
#> bullentin_type: reanalysis
#> cmems_product_id: BALTICSEA_REANALYSIS_PHY_003_011
#> title: CMEMS V4 Reanalysis: NEMO model 3D fields (monthly means)
#> FROM_ORIGINAL_FILE__easternmost_longitude: 30.2357654571533
#> FROM_ORIGINAL_FILE__northernmost_latitude: 65.8914184570312
#> FROM_ORIGINAL_FILE__westernmost_longitude: 9.01375484466553
#> FROM_ORIGINAL_FILE__southernmost_latitude: 48.49169921875
#> shallowest_depth: 1.50136542320251
#> deepest_depth: 711.059204101562
#> source: SMHI reanalysis run NORDIC-NS2_1d_20201201_20201201
#> file_quality_index: 1
#> creation_date: 2021-11-09 UTC
#> bullentin_date: 20201201
#> start_date: 2020-12-01 UTC
#> stop_date: 2020-12-01 UTC
#> start_time: 00:00 UTC
#> stop_time: 00:00 UTC
#> _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
# Get longitude and latitude
lon <- ncvar_get(ncin,"longitude")
nlon <- dim(lon)
head(lon)
#> [1] 9.013755 9.069310 9.124865 9.180420 9.235975 9.291530
lat <- ncvar_get(ncin,"latitude")
nlat <- dim(lat)
head(lat)
#> [1] 48.49170 48.52503 48.55836 48.59170 48.62503 48.65836
# Get time
time <- ncvar_get(ncin,"time")
time
#> [1] 15721.5 15751.0 15780.5 15811.0 15841.5 15872.0 15902.5 15933.5 15964.0
#> [10] 15994.5 16025.0 16055.5 16086.5 16116.0 16145.5 16176.0 16206.5 16237.0
#> [19] 16267.5 16298.5 16329.0 16359.5 16390.0 16420.5 16451.5 16481.0 16510.5
#> [28] 16541.0 16571.5 16602.0 16632.5 16663.5 16694.0 16724.5 16755.0 16785.5
#> [37] 16816.5 16846.5 16876.5 16907.0 16937.5 16968.0 16998.5 17029.5 17060.0
#> [46] 17090.5 17121.0 17151.5 17182.5 17212.0 17241.5 17272.0 17302.5 17333.0
#> [55] 17363.5 17394.5 17425.0 17455.5 17486.0 17516.5 17547.5 17577.0 17606.5
#> [64] 17637.0 17667.5 17698.0 17728.5 17759.5 17790.0 17820.5 17851.0 17881.5
#> [73] 17912.5 17942.0 17971.5 18002.0 18032.5 18063.0 18093.5 18124.5 18155.0
#> [82] 18185.5 18216.0 18246.5 18277.5 18307.5 18337.5 18368.0 18398.5 18429.0
#> [91] 18459.5 18490.5 18521.0 18551.5 18582.0 18612.5 18643.5 18673.0 18702.5
#> [100] 18733.0 18763.5 18794.0 18824.5 18855.5 18886.0 18916.5 18947.0 18977.5
#> [109] 19008.5 19038.0 19067.5 19098.0 19128.5 19159.0 19189.5 19220.5 19251.0
#> [118] 19281.5 19312.0 19342.5 19373.5 19403.0 19432.5 19463.0 19493.5 19524.0
#> [127] 19554.5 19585.5 19616.0 19646.5 19677.0 19707.5 19738.5 19768.5 19798.5
#> [136] 19829.0 19859.5 19890.0 19920.5 19951.5 19982.0 20012.5 20043.0 20073.5
#> [145] 20104.5 20134.0 20163.5 20194.0 20224.5 20255.0 20285.5 20316.5 20347.0
#> [154] 20377.5 20408.0 20438.5 20469.5 20499.0 20528.5 20559.0 20589.5 20620.0
#> [163] 20650.5 20681.5 20712.0 20742.5 20773.0 20803.5 20834.5 20864.0 20893.5
#> [172] 20924.0 20954.5 20985.0 21015.5 21046.5 21077.0 21107.5 21138.0 21168.5
#> [181] 21199.5 21229.5 21259.5 21290.0 21320.5 21351.0 21381.5 21412.5 21443.0
#> [190] 21473.5 21504.0 21534.5 21565.5 21595.0 21624.5 21655.0 21685.5 21716.0
#> [199] 21746.5 21777.5 21808.0 21838.5 21869.0 21899.5 21930.5 21960.0 21989.5
#> [208] 22020.0 22050.5 22081.0 22111.5 22142.5 22173.0 22203.5 22234.0 22264.5
#> [217] 22295.5 22325.0 22354.5 22385.0 22415.5 22446.0 22476.5 22507.5 22538.0
#> [226] 22568.5 22599.0 22629.5 22660.5 22690.5 22720.5 22751.0 22781.5 22812.0
#> [235] 22842.5 22873.5 22904.0 22934.5 22965.0 22995.5 23026.5 23056.0 23085.5
#> [244] 23116.0 23146.5 23177.0 23207.5 23238.5 23269.0 23299.5 23330.0 23360.5
#> [253] 23391.5 23421.0 23450.5 23481.0 23511.5 23542.0 23572.5 23603.5 23634.0
#> [262] 23664.5 23695.0 23725.5 23756.5 23786.0 23815.5 23846.0 23876.5 23907.0
#> [271] 23937.5 23968.5 23999.0 24029.5 24060.0 24090.5 24121.5 24151.5 24181.5
#> [280] 24212.0 24242.5 24273.0 24303.5 24334.5 24365.0 24395.5 24426.0 24456.5
#> [289] 24487.5 24517.0 24546.5 24577.0 24607.5 24638.0 24668.5 24699.5 24730.0
#> [298] 24760.5 24791.0 24821.5 24852.5 24882.0 24911.5 24942.0 24972.5 25003.0
#> [307] 25033.5 25064.5 25095.0 25125.5 25156.0 25186.5 25217.5 25247.0 25276.5
#> [316] 25307.0 25337.5 25368.0 25398.5 25429.5 25460.0 25490.5 25521.0 25551.5
#> [325] 25582.5 25612.5 25642.5 25673.0 25703.5 25734.0 25764.5 25795.5 25826.0
#> [334] 25856.5 25887.0 25917.5
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
#> [1] 336
tunits
#> $hasatt
#> [1] TRUE
#>
#> $value
#> [1] "days since 1950-01-01 00:00:00"
# Get temperature
dname <- "bottomT"
temp_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(temp_array)
#> [1] 383 523 336
# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")
# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])
# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))
# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)
# Replace netCDF fill values with NA's
temp_array[temp_array == fillvalue$value] <- NA
# Next, we need to work with the months that correspond to the quarters that we use.
# loop through each time step, and if it is a good month save it as a raster.
# First get the index of months that correspond to Q4
months
#> [1] 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1
#> [26] 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2
#> [51] 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3
#> [76] 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4
#> [101] 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5
#> [126] 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6
#> [151] 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7
#> [176] 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8
#> [201] 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9
#> [226] 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10
#> [251] 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11
#> [276] 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12
#> [301] 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1
#> [326] 2 3 4 5 6 7 8 9 10 11 12
index_keep_q1 <- which(months < 4)
index_keep_q4 <- which(months > 9)
temp_q1 <- temp_array[, , index_keep_q1]
temp_q4 <- temp_array[, , index_keep_q4]
months_keep_q1 <- months[index_keep_q1]
months_keep_q4 <- months[index_keep_q4]
years_keep_q1 <- years[index_keep_q1]
years_keep_q4 <- years[index_keep_q4]
# Now we have an array with data for that quarter
# We need to now calculate the average within a year.
# Get a sequence that takes every third value between 1: number of months (length)
loop_seq_q1 <- seq(1, dim(temp_q1)[3], by = 3)
loop_seq_q4 <- seq(1, dim(temp_q4)[3], by = 3)
# Create objects that will hold data
dlist_q1 <- list()
dlist_q4 <- list()
temp_1 <- c()
temp_2 <- c()
temp_3 <- c()
temp_ave_q1 <- c()
temp_10 <- c()
temp_11 <- c()
temp_12 <- c()
temp_ave_q4 <- c()
# Now average by quarter. The vector loop_seq_q1 is 1, 4, 7 etc. So first i is 1, 2, 3,
# which is the index we want.
for(i in loop_seq_q1) {
temp_1 <- temp_q1[, , (i)]
temp_2 <- temp_q1[, , (i + 1)]
temp_3 <- temp_q1[, , (i + 2)]
temp_10 <- temp_q4[, , (i)]
temp_11 <- temp_q4[, , (i + 1)]
temp_12 <- temp_q4[, , (i + 2)]
temp_ave_q1 <- (temp_1 + temp_2 + temp_3) / 3
temp_ave_q4 <- (temp_10 + temp_11 + temp_12) / 3
list_pos_q1 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)
list_pos_q4 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)
dlist_q1[[list_pos_q1]] <- temp_ave_q1
dlist_q4[[list_pos_q4]] <- temp_ave_q4
}
# Now name the lists with the year:
names(dlist_q1) <- unique(years_keep_q1)
names(dlist_q4) <- unique(years_keep_q4)
# Now I need to make a loop where I extract the raster value for each year...
# The data is called dat so far in this script
# Filter years in the data frame to only have the years I have temperature for
d_sub_temp_q1 <- dat %>% filter(quarter == 1) %>% filter(year %in% names(dlist_q1)) %>% droplevels()
#> filter: removed 457,436 rows (50%), 457,436 rows remaining
#> filter: no rows removed
d_sub_temp_q4 <- dat %>% filter(quarter == 4) %>% filter(year %in% names(dlist_q4)) %>% droplevels()
#> filter: removed 457,436 rows (50%), 457,436 rows remaining
#> filter: no rows removed
# Create data holding object
temp_data_list_q1 <- list()
temp_data_list_q4 <- list()
# Create factor year for indexing the list in the loop
d_sub_temp_q1$year_f <- as.factor(d_sub_temp_q1$year)
d_sub_temp_q4$year_f <- as.factor(d_sub_temp_q4$year)
# Loop through each year and extract raster values for the data points
for(i in unique(d_sub_temp_q1$year_f)) {
# Set plot limits
ymin = 54; ymax = 59; xmin = 12; xmax = 22
# Subset a year
temp_slice_q1 <- dlist_q1[[i]]
temp_slice_q4 <- dlist_q4[[i]]
# Create raster for that year (i)
r_q1 <- raster(t(temp_slice_q1), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
r_q4 <- raster(t(temp_slice_q4), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
# Flip...
r_q1 <- flip(r_q1, direction = 'y')
r_q4 <- flip(r_q4, direction = 'y')
plot(r_q1, main = paste(i, "Q1"))
plot(r_q4, main = paste(i, "Q4"))
# Change projection to UTM (same as pred grid)
sr <- "+proj=utm +zone=33 +datum=WGS84 +units=m "
proj_raster_q1 <- projectRaster(r_q1, crs = sr)
proj_raster_q4 <- projectRaster(r_q4, crs = sr)
# Filter the same years and select only coordinates
d_slice_q1 <- d_sub_temp_q1 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)
d_slice_q4 <- d_sub_temp_q4 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)
# Make into a SpatialPoints object
data_sp_q1 <- SpatialPoints(d_slice_q1)
data_sp_q4 <- SpatialPoints(d_slice_q4)
# Extract raster value (temperature)
rasValue_q1 <- raster::extract(proj_raster_q1, data_sp_q1, method = "bilinear")
rasValue_q4 <- raster::extract(proj_raster_q4, data_sp_q4, method = "bilinear")
# Now we want to plot the results of the raster extractions by plotting the
# data points over a raster and saving it for each year.
# Make the SpatialPoints object into a raster again (for pl)
df_q1 <- as.data.frame(data_sp_q1)
df_q4 <- as.data.frame(data_sp_q4)
# Add in the raster value in the df holding the coordinates for the data
d_slice_q1$temp <- rasValue_q1
d_slice_q4$temp <- rasValue_q4
# Add in which year
d_slice_q1$year <- i
d_slice_q4$year <- i
# Create a index for the data last where we store all years (because our loop index
# i is not continuous, we can't use it directly)
index_q1 <- as.numeric(d_slice_q1$year)[1] - 1992
index_q4 <- as.numeric(d_slice_q4$year)[1] - 1992
# Add each years' data in the list
temp_data_list_q1[[index_q1]] <- d_slice_q1
temp_data_list_q4[[index_q4]] <- d_slice_q4
}
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
# Now create a data frame from the list of all annual values
big_dat_temp_q1 <- dplyr::bind_rows(temp_data_list_q1)
big_dat_temp_q4 <- dplyr::bind_rows(temp_data_list_q4)
big_dat_temp <- bind_rows(mutate(big_dat_temp_q1, quarter = 1),
mutate(big_dat_temp_q4, quarter = 4))
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA
#> mutate: new variable 'quarter' (double) with one unique value and 0% NASalinity
# https://data.marine.copernicus.eu/product/BALTICSEA_REANALYSIS_PHY_003_011/download?dataset=dataset-reanalysis-nemo-monthlymeans
# Open the netCDF file
ncin <- nc_open("data/NEMO_Nordic_SCOBI/dataset-reanalysis-nemo-monthlymeans_1668587452211.nc")
print(ncin)
#> File data/NEMO_Nordic_SCOBI/dataset-reanalysis-nemo-monthlymeans_1668587452211.nc (NC_FORMAT_CLASSIC):
#>
#> 1 variables (excluding dimension variables):
#> float sob[longitude,latitude,time]
#> long_name: Sea water salinity at sea floor
#> missing_value: NaN
#> standard_name: sea_water_salinity
#> units: 0.001
#> _FillValue: NaN
#> _ChunkSizes: 1
#> _ChunkSizes: 523
#> _ChunkSizes: 383
#>
#> 3 dimensions:
#> time Size:335
#> axis: T
#> long_name: Validity time
#> standard_name: time
#> units: days since 1950-01-01 00:00:00
#> calendar: gregorian
#> _ChunkSizes: 512
#> _CoordinateAxisType: Time
#> valid_min: 15751
#> valid_max: 25917.5
#> latitude Size:187
#> axis: Y
#> standard_name: latitude
#> long_name: latitude
#> units: degrees_north
#> _CoordinateAxisType: Lat
#> valid_min: 53.1249580383301
#> valid_max: 59.3248596191406
#> longitude Size:199
#> standard_name: longitude
#> long_name: longitude
#> units: degrees_east
#> axis: X
#> _CoordinateAxisType: Lon
#> valid_min: 11.1248445510864
#> valid_max: 22.12473487854
#>
#> 24 global attributes:
#> references: http://www.smhi.se
#> institution: Swedish Meterological and Hydrological Institute
#> history: See source and creation_date attributees
#> Conventions: CF-1.5
#> contact: servicedesk_cmems@mercator-ocean.eu
#> comment: Provided by SMHI as a Copernicus Marine Environment Monitoring Service production unit
#> bullentin_type: reanalysis
#> cmems_product_id: BALTICSEA_REANALYSIS_PHY_003_011
#> title: CMEMS V4 Reanalysis: NEMO model 3D fields (monthly means)
#> FROM_ORIGINAL_FILE__easternmost_longitude: 30.2357654571533
#> FROM_ORIGINAL_FILE__northernmost_latitude: 65.8914184570312
#> FROM_ORIGINAL_FILE__westernmost_longitude: 9.01375484466553
#> FROM_ORIGINAL_FILE__southernmost_latitude: 48.49169921875
#> shallowest_depth: 1.50136542320251
#> deepest_depth: 711.059204101562
#> source: SMHI reanalysis run NORDIC-NS2_1d_20201201_20201201
#> file_quality_index: 1
#> creation_date: 2021-11-09 UTC
#> bullentin_date: 20201201
#> start_date: 2020-12-01 UTC
#> stop_date: 2020-12-01 UTC
#> start_time: 00:00 UTC
#> stop_time: 00:00 UTC
#> _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
# Get longitude and latitude
lon <- ncvar_get(ncin,"longitude")
nlon <- dim(lon)
head(lon)
#> [1] 11.12484 11.18040 11.23596 11.29151 11.34706 11.40262
lat <- ncvar_get(ncin,"latitude")
nlat <- dim(lat)
head(lat)
#> [1] 53.12496 53.15829 53.19162 53.22496 53.25829 53.29162
# Get time
time <- ncvar_get(ncin,"time")
time
#> [1] 15751.0 15780.5 15811.0 15841.5 15872.0 15902.5 15933.5 15964.0 15994.5
#> [10] 16025.0 16055.5 16086.5 16116.0 16145.5 16176.0 16206.5 16237.0 16267.5
#> [19] 16298.5 16329.0 16359.5 16390.0 16420.5 16451.5 16481.0 16510.5 16541.0
#> [28] 16571.5 16602.0 16632.5 16663.5 16694.0 16724.5 16755.0 16785.5 16816.5
#> [37] 16846.5 16876.5 16907.0 16937.5 16968.0 16998.5 17029.5 17060.0 17090.5
#> [46] 17121.0 17151.5 17182.5 17212.0 17241.5 17272.0 17302.5 17333.0 17363.5
#> [55] 17394.5 17425.0 17455.5 17486.0 17516.5 17547.5 17577.0 17606.5 17637.0
#> [64] 17667.5 17698.0 17728.5 17759.5 17790.0 17820.5 17851.0 17881.5 17912.5
#> [73] 17942.0 17971.5 18002.0 18032.5 18063.0 18093.5 18124.5 18155.0 18185.5
#> [82] 18216.0 18246.5 18277.5 18307.5 18337.5 18368.0 18398.5 18429.0 18459.5
#> [91] 18490.5 18521.0 18551.5 18582.0 18612.5 18643.5 18673.0 18702.5 18733.0
#> [100] 18763.5 18794.0 18824.5 18855.5 18886.0 18916.5 18947.0 18977.5 19008.5
#> [109] 19038.0 19067.5 19098.0 19128.5 19159.0 19189.5 19220.5 19251.0 19281.5
#> [118] 19312.0 19342.5 19373.5 19403.0 19432.5 19463.0 19493.5 19524.0 19554.5
#> [127] 19585.5 19616.0 19646.5 19677.0 19707.5 19738.5 19768.5 19798.5 19829.0
#> [136] 19859.5 19890.0 19920.5 19951.5 19982.0 20012.5 20043.0 20073.5 20104.5
#> [145] 20134.0 20163.5 20194.0 20224.5 20255.0 20285.5 20316.5 20347.0 20377.5
#> [154] 20408.0 20438.5 20469.5 20499.0 20528.5 20559.0 20589.5 20620.0 20650.5
#> [163] 20681.5 20712.0 20742.5 20773.0 20803.5 20834.5 20864.0 20893.5 20924.0
#> [172] 20954.5 20985.0 21015.5 21046.5 21077.0 21107.5 21138.0 21168.5 21199.5
#> [181] 21229.5 21259.5 21290.0 21320.5 21351.0 21381.5 21412.5 21443.0 21473.5
#> [190] 21504.0 21534.5 21565.5 21595.0 21624.5 21655.0 21685.5 21716.0 21746.5
#> [199] 21777.5 21808.0 21838.5 21869.0 21899.5 21930.5 21960.0 21989.5 22020.0
#> [208] 22050.5 22081.0 22111.5 22142.5 22173.0 22203.5 22234.0 22264.5 22295.5
#> [217] 22325.0 22354.5 22385.0 22415.5 22446.0 22476.5 22507.5 22538.0 22568.5
#> [226] 22599.0 22629.5 22660.5 22690.5 22720.5 22751.0 22781.5 22812.0 22842.5
#> [235] 22873.5 22904.0 22934.5 22965.0 22995.5 23026.5 23056.0 23085.5 23116.0
#> [244] 23146.5 23177.0 23207.5 23238.5 23269.0 23299.5 23330.0 23360.5 23391.5
#> [253] 23421.0 23450.5 23481.0 23511.5 23542.0 23572.5 23603.5 23634.0 23664.5
#> [262] 23695.0 23725.5 23756.5 23786.0 23815.5 23846.0 23876.5 23907.0 23937.5
#> [271] 23968.5 23999.0 24029.5 24060.0 24090.5 24121.5 24151.5 24181.5 24212.0
#> [280] 24242.5 24273.0 24303.5 24334.5 24365.0 24395.5 24426.0 24456.5 24487.5
#> [289] 24517.0 24546.5 24577.0 24607.5 24638.0 24668.5 24699.5 24730.0 24760.5
#> [298] 24791.0 24821.5 24852.5 24882.0 24911.5 24942.0 24972.5 25003.0 25033.5
#> [307] 25064.5 25095.0 25125.5 25156.0 25186.5 25217.5 25247.0 25276.5 25307.0
#> [316] 25337.5 25368.0 25398.5 25429.5 25460.0 25490.5 25521.0 25551.5 25582.5
#> [325] 25612.5 25642.5 25673.0 25703.5 25734.0 25764.5 25795.5 25826.0 25856.5
#> [334] 25887.0 25917.5
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
#> [1] 335
tunits
#> $hasatt
#> [1] TRUE
#>
#> $value
#> [1] "days since 1950-01-01 00:00:00"
# Get Salinity
dname <- "sob"
sal_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(sal_array)
#> [1] 199 187 335
# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")
# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])
# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))
# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)
# Replace netCDF fill values with NA's
sal_array[sal_array == fillvalue$value] <- NA
# Next, we need to work with the months that correspond to the quarters that we use.
# loop through each time step, and if it is a good month save it as a raster.
# First get the index of months that correspond to Q4
months
#> [1] 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2
#> [26] 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3
#> [51] 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4
#> [76] 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5
#> [101] 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6
#> [126] 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7
#> [151] 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8
#> [176] 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9
#> [201] 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10
#> [226] 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11
#> [251] 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12
#> [276] 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1
#> [301] 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2
#> [326] 3 4 5 6 7 8 9 10 11 12
index_keep_q1 <- which(months < 4)
index_keep_q4 <- which(months > 9)
sal_q1 <- sal_array[, , index_keep_q1]
sal_q4 <- sal_array[, , index_keep_q4]
months_keep_q1 <- months[index_keep_q1]
months_keep_q4 <- months[index_keep_q4]
years_keep_q1 <- years[index_keep_q1]
years_keep_q4 <- years[index_keep_q4]
# Now we have an array with data for that quarter
# We need to now calculate the average within a year.
# Get a sequence that takes every third value between 1: number of months (length)
loop_seq_q1 <- seq(1, dim(sal_q1)[3], by = 3)
loop_seq_q4 <- seq(1, dim(sal_q4)[3], by = 3)
# Create objects that will hold data
dlist_q1 <- list()
dlist_q4 <- list()
sal_1 <- c()
sal_2 <- c()
sal_3 <- c()
sal_ave_q1 <- c()
sal_10 <- c()
sal_11 <- c()
sal_12 <- c()
sal_ave_q4 <- c()
# Now average by quarter. The vector loop_seq_q1 is 1, 4, 7 etc. So first i is 1, 2, 3,
# which is the index we want.
dim(sal_q1)
#> [1] 199 187 83
dim(sal_q4)
#> [1] 199 187 84
# Hmm, we didn't get the first month in the salinity series... repeat month 2 and fill in so the dimensions are correct
sal_q1 <- sal_q1[,,c(1, 1:83)]
dim(sal_q1)
#> [1] 199 187 84
for(i in loop_seq_q1) {
sal_1 <- sal_q1[, , (i)]
sal_2 <- sal_q1[, , (i + 1)]
sal_3 <- sal_q1[, , (i + 2)]
sal_10 <- sal_q4[, , (i)]
sal_11 <- sal_q4[, , (i + 1)]
sal_12 <- sal_q4[, , (i + 2)]
sal_ave_q1 <- (sal_1 + sal_2 + sal_3) / 3
sal_ave_q4 <- (sal_10 + sal_11 + sal_12) / 3
list_pos_q1 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)
list_pos_q4 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)
dlist_q1[[list_pos_q1]] <- sal_ave_q1
dlist_q4[[list_pos_q4]] <- sal_ave_q4
}
# Now name the lists with the year:
names(dlist_q1) <- unique(years_keep_q1)
names(dlist_q4) <- unique(years_keep_q4)
# Now I need to make a loop where I extract the raster value for each year...
# The cpue data is called dat so far in this script
# Filter years in the cpue data frame to only have the years I have salinity for
d_sub_sal_q1 <- dat %>% filter(quarter == 1) %>% filter(year %in% names(dlist_q1)) %>% droplevels()
#> filter: removed 457,436 rows (50%), 457,436 rows remaining
#> filter: no rows removed
d_sub_sal_q4 <- dat %>% filter(quarter == 4) %>% filter(year %in% names(dlist_q4)) %>% droplevels()
#> filter: removed 457,436 rows (50%), 457,436 rows remaining
#> filter: no rows removed
# Create data holding object
sal_data_list_q1 <- list()
sal_data_list_q4 <- list()
# Create factor year for indexing the list in the loop
d_sub_sal_q1$year_f <- as.factor(d_sub_sal_q1$year)
d_sub_sal_q4$year_f <- as.factor(d_sub_sal_q4$year)
# Loop through each year and extract raster values for the cpue data points
for(i in unique(d_sub_sal_q1$year_f)) {
# Set plot limits
ymin = 54; ymax = 58; xmin = 12; xmax = 22
# Subset a year
sal_slice_q1 <- dlist_q1[[i]]
sal_slice_q4 <- dlist_q4[[i]]
# Create raster for that year (i)
r_q1 <- raster(t(sal_slice_q1), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
r_q4 <- raster(t(sal_slice_q4), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
# Flip...
r_q1 <- flip(r_q1, direction = 'y')
r_q4 <- flip(r_q4, direction = 'y')
plot(r_q1, main = paste(i, "Q1"))
plot(r_q4, main = paste(i, "Q4"))
# Change projection to UTM (same as pred grid)
sr <- "+proj=utm +zone=33 +datum=WGS84 +units=m "
proj_raster_q1 <- projectRaster(r_q1, crs = sr)
proj_raster_q4 <- projectRaster(r_q4, crs = sr)
# Filter the same year (i) in the cpue data and select only coordinates
d_slice_q1 <- d_sub_sal_q1 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)
d_slice_q4 <- d_sub_sal_q4 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)
# Make into a SpatialPoints object
data_sp_q1 <- SpatialPoints(d_slice_q1)
data_sp_q4 <- SpatialPoints(d_slice_q4)
# Extract raster value (salinity)
rasValue_q1 <- raster::extract(proj_raster_q1, data_sp_q1, method = "bilinear")
rasValue_q4 <- raster::extract(proj_raster_q4, data_sp_q4, method = "bilinear")
# Now we want to plot the results of the raster extractions by plotting the cpue
# data points over a raster and saving it for each year.
# Make the SpatialPoints object into a raster again (for pl)
df_q1 <- as.data.frame(data_sp_q1)
df_q4 <- as.data.frame(data_sp_q4)
# Add in the raster value in the df holding the coordinates for the cpue data
d_slice_q1$sal <- rasValue_q1
d_slice_q4$sal <- rasValue_q4
# Add in which year
d_slice_q1$year <- i
d_slice_q4$year <- i
# Create a index for the data last where we store all years (because our loop index
# i is not continuous, we can't use it directly)
index_q1 <- as.numeric(d_slice_q1$year)[1] - 1992
index_q4 <- as.numeric(d_slice_q4$year)[1] - 1992
# Add each years' data in the list
sal_data_list_q1[[index_q1]] <- d_slice_q1
sal_data_list_q4[[index_q4]] <- d_slice_q4
}
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#> changed 16,337 values (100%) of 'Y' (0 new NA)
# Now create a data frame from the list of all annual values
big_dat_sal_q1 <- dplyr::bind_rows(sal_data_list_q1)
big_dat_sal_q4 <- dplyr::bind_rows(sal_data_list_q4)
big_dat_sal <- bind_rows(mutate(big_dat_sal_q1, quarter = 1),
mutate(big_dat_sal_q4, quarter = 4))
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA
#> mutate: new variable 'quarter' (double) with one unique value and 0% NAbig_dat_oxy <- big_dat_oxy %>% mutate(id_env = paste(year, quarter, X, Y, sep = "_"))
#> mutate: new variable 'id_env' (character) with 914,872 unique values and 0% NA
big_dat_temp <- big_dat_temp %>% mutate(id_env = paste(year, quarter, X, Y, sep = "_"))
#> mutate: new variable 'id_env' (character) with 914,872 unique values and 0% NA
big_dat_sal <- big_dat_sal %>% mutate(id_env = paste(year, quarter, X, Y, sep = "_"))
#> mutate: new variable 'id_env' (character) with 914,872 unique values and 0% NA
env_dat <- left_join(big_dat_oxy, big_dat_temp) %>%
left_join(big_dat_sal) %>%
dplyr::select(id_env, oxy, temp, sal)
#> Joining, by = c("X", "Y", "year", "quarter", "id_env")
#> left_join: added one column (temp)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 914,872
#> > =========
#> > rows total 914,872
#> Joining, by = c("X", "Y", "year", "quarter", "id_env")
#> left_join: added one column (sal)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 914,872
#> > =========
#> > rows total 914,872
dat <- dat %>%
mutate(id_env = paste(year, quarter, X*1000, Y*1000, sep = "_")) %>%
left_join(env_dat) %>%
dplyr::select(-id_env)
#> mutate: new variable 'id_env' (character) with 914,872 unique values and 0% NA
#> Joining, by = "id_env"left_join: added 3 columns (oxy, temp, sal)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 914,872
#> > =========
#> > rows total 914,872
dat <- dat %>% drop_na(temp) %>% drop_na(oxy) %>% drop_na(sal)
#> drop_na: removed 5,712 rows (1%), 909,160 rows remaining
#> drop_na: no rows removed
#> drop_na: removed 1,008 rows (<1%), 908,152 rows remaining
# Temperature
plot_map_fc +
geom_raster(data = filter(dat, quarter == 4), aes(X*1000, Y*1000, fill = temp)) +
facet_wrap(~year)
#> filter: removed 454,076 rows (50%), 454,076 rows remaining
#> Warning: Removed 21336 rows containing missing values (geom_raster).
# Oxygen
plot_map_fc +
geom_raster(data = filter(dat, quarter == 4), aes(X*1000, Y*1000, fill = oxy)) +
facet_wrap(~year)
#> filter: removed 454,076 rows (50%), 454,076 rows remaining
#> Warning: Removed 21336 rows containing missing values (geom_raster).
# Salinity
plot_map_fc +
geom_raster(data = filter(dat, quarter == 4), aes(X*1000, Y*1000, fill = sal)) +
facet_wrap(~year)
#> filter: removed 454,076 rows (50%), 454,076 rows remaining
#> Warning: Removed 21336 rows containing missing values (geom_raster).
pred_grid <- datAdd ICES areas
# https://stackoverflow.com/questions/34272309/extract-shapefile-value-to-point-with-r
# https://gis.ices.dk/sf/
shape <- shapefile("data/ICES_StatRec_mapto_ICES_Areas/StatRec_map_Areas_Full_20170124.shp")
head(shape)
#> OBJECTID ID ICESNAME SOUTH WEST NORTH EAST AREA_KM2 stat_x stat_y Area_27
#> 0 1 47 47A0 59.0 -44 59.5 -43 3178 -43.5 59.25 14.b.2
#> 1 2 48 48A0 59.5 -44 60.0 -43 3132 -43.5 59.75 14.b.2
#> 2 3 49 49A0 60.0 -44 60.5 -43 3085 -43.5 60.25 14.b.2
#> 3 4 50 50A0 60.5 -44 61.0 -43 3038 -43.5 60.75 14.b.2
#> 4 5 51 51A0 61.0 -44 61.5 -43 2991 -43.5 61.25 14.b.2
#> 5 6 52 52A0 61.5 -44 62.0 -43 2944 -43.5 61.75 14.b.2
#> Perc MaxPer RNDMaxPer AreasList Shape_Leng Shape_Area
#> 0 100.00000000 100.0000000 100 14.b.2 3 0.5
#> 1 84.12674145 84.1267414 84 14.b.2 3 0.5
#> 2 24.99803694 24.9980369 25 14.b.2 3 0.5
#> 3 11.97744244 11.9774424 12 14.b.2 3 0.5
#> 4 3.89717625 3.8971762 4 14.b.2 3 0.5
#> 5 0.28578428 0.2857843 0 14.b.2 3 0.5
pts <- SpatialPoints(cbind(pred_grid$lon, pred_grid$lat),
proj4string = CRS(proj4string(shape)))
#> Warning in proj4string(shape): CRS object has comment, which is lost in output
pred_grid$subdiv <- over(pts, shape)$Area_27
# Rename subdivisions to the more common names and do some more filtering (by sub div and area)
sort(unique(pred_grid$subdiv))
#> [1] "3.d.24" "3.d.25" "3.d.26" "3.d.27" "3.d.28.2"
pred_grid <- pred_grid %>%
mutate(sub_div = factor(subdiv),
sub_div = fct_recode(subdiv,
"24" = "3.d.24",
"25" = "3.d.25",
"26" = "3.d.26",
"27" = "3.d.27",
"28" = "3.d.28.1",
"28" = "3.d.28.2",
"29" = "3.d.29"),
sub_div = as.character(sub_div)) %>%
filter(sub_div %in% c("24", "25", "26", "27", "28", 2)) %>%
filter(lat > 54 & lat < 59 & lon < 22)
#> Warning: Unknown levels in `f`: 3.d.28.1, 3.d.29
#> mutate: new variable 'sub_div' (character) with 5 unique values and 0% NA
#> filter: no rows removed
#> filter: no rows removed
# Add ICES rectangles
pred_grid$ices_rect <- mapplots::ices.rect2(lon = pred_grid$lon, lat = pred_grid$lat)
plot_map +
geom_raster(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000, fill = oxy)) +
facet_wrap(~sub_div)
#> filter: removed 875,718 rows (96%), 32,434 rows remaining
#> Warning: Removed 1524 rows containing missing values (geom_raster).
pred_grid <- pred_grid %>% dplyr::select(-subdiv)Save
# Remove variables and save
pred_grid_93_06 <- pred_grid %>% filter(year < 2007)
#> filter: removed 454,076 rows (50%), 454,076 rows remaining
pred_grid_07_19 <- pred_grid %>% filter(year > 2006)
#> filter: removed 454,076 rows (50%), 454,076 rows remaining
write.csv(pred_grid_93_06, file = "data/clean/pred_grid_(1_2).csv", row.names = FALSE)
write.csv(pred_grid_07_19, file = "data/clean/pred_grid_(2_2).csv", row.names = FALSE)